Grain Boundary Diffusion in Copper under Tensile Stress 
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Stress enhanced self-diffusion of Copper on the E3 twin grain boundary was examined with molec- 
ular dynamics simulations. The presence of uniaxial tensile stress results in a significant reduction 
in activation energy for grain-boundary self-diffusion of magnitude 5 eV per unit strain. Using a 
theoretical model of point defect formation and diffusion, the functional dependence of the effective 
activation energy Q on uniaxial tensile strain e is shown to be described by Q{e) — Qo — EoV*e 
where Eo is the zero-temperature Young's modulus and V* is an effective activation volume. The 
simulation data agree well with this model and comparison between data and model suggests that 
V* = 0.6Q where fl is the atomic volume. V* /Q, = 0.6 is consistent with a vacancy-dominated 
diffusion mechanism. 



I. INTRODUCTION 

Studies of failure mechanisms in nano-structured 
metallic interconnects show that diffusion near grain 
boundaries (GBs) plays a critical role in initiation of 
destructive mass-transport iSi^ For example, electromi- 
gration failure of poly-crystalline aluminum and copper 
interconnects is initiated predominantly at grain bound- 
aries and triple points. In accord with these phenomeno- 
logical observations are measured and calculated acti- 
vation energies for self-diffusion that are approximately 
one-third to one-half the corresponding bulk activa- 
tion energies.^ Typically metallic interconnects are under 
large tensile strains due to both the deposition process 
and the interaction with an underlying substrate or over- 
lying passivation layer. Metallic interconnects in very- 
and ultra-large-scale-integrated circuit (VLSI/ULSI) ap- 
plications can experience tensile stresses in excess of sev- 
eral hundred MPa resulting in strains on the order of 
several hundredths of a percent. While it is well known 
that external tensile strain enhances diffusivity in bulk 
materials, there have been relatively few experimental or 
computational investigations of diffusion in strained met- 
als despite the technological significance of the problem. 

Due to its use in nano-scale interconnect materi- 
als, the structural properties and self-diffusion mech- 
anisms in copper near GBs have received consider- 
able attentioni^*SiLSi2ii£ Several experiments on the mi- 
crostructure of copper thin films suggest that the S3 
coherent twin grain boundary accounts for roughly 
42% of all coincident-site-lattice (CSL) boundaries in 
as-deposited copper thin filmsiii CSL boundaries to- 
gether represent the majority of GB geometries in poly- 
crystalline copper. Further, calculations of GB surface 
energies show that the E3 boundary has a markedly lower 
energy than all but the E7 boundaryii^ While, compared 
to other CSL boundaries, defect mobility is relatively 
low in the E3 GB, the large number of such boundaries 
present in copper thin films motivates the study of diffu- 
sion near these structures. 

Bulk self-diffusion in metals may be considerably sim- 
pler than diffusion near GBs. Sandberg et al. have re- 
cently demonstrated that monovacancies are the only ac- 



tive contributors to bulk diffusion in aluminum at all tem- 
peratures up to the melting point. The situation near 
GBs is far from resolved, however, and is certainly more 
complicated than that in the bulk. Sorensen et al. es- 
tablished that symmetric tilt GBs support both vacancy 
and interstitial-mediated diffusion at low temperatures, 
while at elevated temperatures, more complex defect in- 
teractions such as Frenkel pairs and even ring mecha- 
nisms contribute to GB diffusion.^ Recently, Suzuki and 
Mishin demonstrated that vacancy and interstitial mech- 
anisms can have comparable formation energies in sym- 
metric tilt boundaries in copperii^ There is also ample 
experimental evidence of complex, cooperative motion of 
large numbers of atoms in both twist and tilt GBs in thin 
Au films at high temperatureJ^ 

Nomura and Adams have extensively studied diffu- 
sion mechanisms in and near both symmetric tilt and 
pure twist GBs in copper i4 Assuming only a vacancy 
mechanism for diffusion in twist boundaries, Nomura and 
Adams found that diffusion at low temperatures is dom- 
inated by migration along the screw dislocations com- 
prising high angle twist boundaries, while high tempera- 
ture diffusion occurred primarily through the bulk. No- 
mura and others have also established that vacancies are 
strongly bound to the GB itself so that vacancy diffusion 
is primarily restricted to the plane defined by the GB. 

Because of the constraining geometry of the twin 
boundary studied in the present work, it is unlikely that 
interstitials and more complex defects such as Frenkel 
pairs contribute significantly to diffusion on the S3 GB 
at low temperatures!^ At elevated temperatures, how- 
ever, defect concentrations may be high enough that in- 
teractions between point defects may contribute to dif- 
fusion near the GB. While appropriate mechanisms have 
been tentatively identified for bulk diffusion, mechanisms 
for diffusion near grain boundaries have not been clearly 
established 

As a computational problem, diffusion near GBs in 
copper has been examined in the context of molecular 
dynamics (MD), molecular statics (MS), and, more re- 
cently, kinetic Monte Carlo (KMC) models.^'^'^'^ Because 
MD evolves the atomic lattice directly according to a 
specified potential, all migration pathways and diffusion 
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mechanisms are allowed. For this reason, MD is useful 
in identifying the kinetics of defect formation, differenti- 
ating between diffusion through vacancy and interstitial 
migration when, as is the case for GB diffusion, a domi- 
nant mechanism has not been identified. 

Typically, MD simulations are limited to time scales 
on the order of nano-seconds, often too short to reliably 
extrapolate diffusion constants at the relatively low (in 
terms of bulk melting temperatures) operating temper- 
atures of copper interconnects in micro-electronics tech- 
nologies. As a result, an accurate estimation of activation 
energies for self-diffusion is difficult. Conversely, the MS 
and KMC methods and other approaches involving the 
coupled application of Monte Carlo and molecular stat- 
ics approaches requires an a priori knowledge of defect 
types, concentrations, and jump frequencies, but over- 
comes the short-time scale limitations of direct MD sim- 
ulations. Recently, hybrid approaches have been demon- 
strated that combine the strengths of each simulation 
technique to more reliably estimate diffusion properties 
in metalsi^ 

In this paper, I report on the use of MD to determine 
activation energies as a function of static lattice-strain for 
self-diffusion on the S3 coherent twin grain boundary in 
copper. In Sec. II, the details of the computational model 
are developed, and fundamental elastic properties of the 
simulated crystal are calculated and compared with exist- 
ing data to validate the bulk mechanical behavior of the 
computational model. A theoretical framework for the 
diffusion data obtained in the simulations is presented in 
Sec. Ill, while in Sec. IV, the activation energy data 
and other relevant results are presented and interpreted 
within the theoretical framework of Sec. III. Finally, a 
discussion of the results in the context of existing and 
future experimental and numerical work is presented in 
Sec. V. 



II. COMPUTATIONAL MODEL 

Molecular dynamics simulations using potentials de- 
rived from the embedded-atom model (EAM) are well- 
known to accurately account for the behaviors of metals, 
particularly near surfaces and point defects where other 
semi-classical potentials fail. The accuracy of several of 
these potentials in accounting for the thermodynamics 
and energetics of extended defects such as grain bound- 
aries has recently been established as welliii The sim- 
ulations reported here make use of EAM potentials for 
copper published by Johnson.^® The potential has been 
modified for a smooth cut-off and to better match bulk 
mechanical properties of copper. The MD time-step is 
0.9 fs, and each simulation is carried out for at least 1.2 
ns, with lower temperature simulations run for up to 2.2 
ns to obtain reasonable diffusion statistics. 

A simulation cell with dimensions 20Ax40Ax40A, 
consisting of roughly 128,000 atoms is sliced in half such 
that the slice-plane separates two {111} crystal faces 




FIG. 1: A 12Ax25Ax25A section of the full simulation cell 
with a E3 GB on the central (111) plane. A tensile strain 
is applied along the external coordinate direction x which is 
aligned along the (111) crystallographic direction. 



of the copper lattice. The slice plane defines a grain 
boundary of area 1600A^. To generate the coherent twin 
boundary, one side of the crystal is rotated by 60° about 
the (111) axis. The geometry is illustrated in Fig. ^ Pe- 
riodic boundary conditions are applied to the faces of the 
crystal normal to the GB, while the two surfaces parallel 
to the GB are free. As indicated in Fig. ^ the exter- 
nal coordinate direction x is aligned with the crystallo- 
graphic direction (111). 

To find the energetically favorable separation of the 
two crystallites defining the GB, Monte Carlo relaxation 
is carried out to optimize the separation between the two 
crystallites and a minimum-energy separation is obtained 
for each temperature. The equilibrium lattice constant is 
established through short MD calculations which search 
for a minimum in internal stress as the lattice constant 
is varied at a fixed temperature. 

For temperatures in the range OAATm <T < 0.7ATm, 
where T,„ — 1356K is the experimental bulk melting tem- 
perature of copper, the Young's modulus, E = au/cu 
and the Poisson ratio, v = £22/611 are evaluated from the 
elastic coefficients Cn and C12. The Cij are obtained 
from the ensemble averages of fluctuations in the stress 
tensor using the method of Ref. 20. Poisson's ratio v 
for tension in the x direction is relatively insensitive to 
temperature and takes the value v = 0.31 ± 0.02 for all 
simulations reported here. 

The calculated temperature dependence of the Young's 
modulus E{T) is shown in Fig. |21 and obeys the phe- 
nomenological relation 

E{T)^Eoil-jT/Tra) (1) 

for temperatures T < O.TiTm where Eq = 116 GPa is 
the zero-temperature modulus, and 7 « 0.55. Qagin 
et al. calculate the elastic constant Cn as a function 
of temperature for copper. Their data is in excellent 
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FIG. 2: The temperature dependence of the Young's modulus 
E = ail /en for uniaxial tension in units of the zero temper- 
ature modulus Eo — 116 GPa (circles). The solid line is a 
linear regression fit with slope —0.55. The triangles are data 
obtained from Ref,^- with Eq — 127 GPa. Both data sets are 
well described by E{T) = Eo{l - -yT/Tm) with 7 = 0.55 and 
Tm = 1356ii:. 



agreement with Eq. ^ when one applies the relation 
Cii = (1 - i^)E/{l + z/)(l - 2v) to obtain the Young's 
modulus for uniaxial tension along the x direction from 
the data in Refi^. 

To study diffusion under strain, the lattice is equili- 
brated at each target temperature using Monte Carlo re- 
laxation. A strain rate of e = 10^'*ps~^ is applied to the 
(111) faces of the lattice until a desired average internal 
static strain en is reached. Diffusivities are calculated 
from the mean-square displacement of atoms in and near 
the GB according to 



(2) 



where ra{t) is the a component of the atomic displace- 
ment vector on the GB plane at time t, and the average 
is over all atoms originating in either of the two atomic 
planes defining the grain boundary. The diffusivity D|| 
parallel to the GB is computed at temperatures ranging 
from OAATm, to 0.74Tm for a range of static strains from 
to 0.07. The maximum corresponding stresses are on 
the order of a few GPa. From the diffusivity data, Ar- 
rhenius plots are constructed and the effective activation 
energy is computed at each value of the static strain e. 



III. DIFFUSION UNDER UNIAXIAL TENSILE 
STRAIN 

In the simulations reported here, the copper crystal is 
under uniaxial tension an due to a static strain applied 
at the (111) faces as shown in Fig. ^ The lattice strain is 



given by en = e, €22 = £33 = —ve, and all other e^ 
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diffusivity components D12 and D13 which we can take 
to be equal: Z^n = D12 = -D13. Transport normal to 
the GB is described by the diffusion constant Dn. In 
the simulations reported in this paper, we have adopted 
a very narrow definition of the grain boundary region 
as consisting of just the two atomic planes adjacent to 
the interface. For this reason, we are interested in the 
diffusivity D\\ within the grain boundary plane y — z, 
and we anticipate that Dn will reflect bulk diffusivity 
that is not of interest here. 

Experimentally, the effect of tensile strain e on diffu- 
sivity D has been characterized by an effective reduction 
of the activation energy for diffusion by the amount^SiS^ 



g' = -kT 



d\nD 



(3) 



Aziz has provided a rigorous theoretical foundation for 
the empirical result of Eq. l^jSiS^ We make use of the 
theory in Reim^ in what follows. 

After Aziz, we consider the components of the diffu- 
sivity tensor D to be proportional to the product of the 
defect concentration C (cr) and appropriate component of 
the mobility tensor M(cr) of point defects in the lattice. 
Here, cr is the spatially uniform stress tensor with com- 
ponents aij . The equilibrium concentration of defects is 
related to the stress- free concentration C(0) by 



Cja) _ exp (cr ■ yf) 

C(0) ~ kT 



(4) 



The formation strain tensor V-^ depends on the defect 
mechanism, as well as on the location of the defect source 
and represents the volume changes due to the formation 
of a defect. Experimentally, only certain components of 
are accessible. In the case of hydrostatic stresses, the 
formation volume V-^ = TrV^ is measured. 

Consider point defects forming on or near the GB 
plane. A vacancy (interstitial) defect initially increases 
(decreases) the sample volume by one atomic volume unit 
n on the (111) face which is normal to the external coor- 
dinate direction x. Subsequent to a vacancy (interstitial) 
formation, the lattice relaxes around the defect, contract- 
ing (expanding) around the vacancy (interstitial). The 
resulting change in volume is characterized by the relax- 
ation volume which is negative (positive) for vacancy 
(interstitial) defects. The defect formation tensor is then 
given by 




(5) 



Diffusion in the plane of the GB is characterized by the 



where + (— ) sign is for vacancy (interstitial) formation. 

The directional dependence of the diffusivity is con- 
tained in the mobility tensor, M. The defect mobility 
is also stress dependent, and its contribution to the self- 
diffusivity is through the migration strain tensor V™. 
Diffusion in the plane of the grain boundary depends on 



the mobility Mii = M12 ~ M13 according to 



M||H 
M||(0) 



cxp 



kT 



(6) 



where V|j" is the migration strain tensor in the plane of 
the GB and the migration volume V™ = TrV™ repre- 
sents the change in volume of the lattice after the defect 
has reached the saddle point in its migration path. 

From Equations 01 and |B1 the stress dependence of _D|| 
is given by 



In- 



Dn 



"^ll(O) 



kT 



(7) 



For the simple tension condition used in our simulations, 
the stress tensor has a single non-zero component, cth = 
Ee where e is the magnitude of the strain along the x 
direction: en — e > 0. Under uniaxial tension in the x 
direction, Eq. [T] becomes 



In 



L'll(e) _ £^(±17 + ^73- 



* XX / 



kT 



(8) 



where is the dimension change in the x direction 
upon defect migration to its saddle point. 

In an unstrained lattice, the usual Arrhenius form for 
the diffusivity, 



Do 



Qo 
kT 



(9) 



holds, where Dq is a material constant, and Qo is the 
activation energy for self-diffusion in the unstrained lat- 
tice. Qo is the sum of a defect formation energy Qq and 
a migration energy, Q™. It is clear from Eqns. [HlandO 
that the effective activation energy for the case of simple 
tension considered here is 



where 



Q(e) = Qo - EV*e, 



(10) 



(11) 



The effect of tensile strain is to reduce the effective ac- 
tivation energy for self-diffusion by the amount EV* per 
unit strain. 

High-temperature elastic softening of the copper crys- 
tal does not change the fundamental Arrhenius analysis 
of diffusion coefhcients. Combining Eqns. Q 13 and 1111 
we find the Diffusion coefficient is still of Arrhenius form: 



In^" 



^llo(e) 



Q(f) 
kT 



(12) 



where D\\ „(e) = D\\oe-^^°^'^/'''^'^ , and Duo is a ma- 
terial constant. The strain-dependent activation energy 
Q(e) in Eq. ^is given by Eq. ^with E replaced by 
Eq. The functional form of the effective activation energy 
is unchanged by elastic softening, and the usual Arrhe- 
nius analysis is still appropriate for the determination of 
diffusivities and activation energies. 




FIG. 3: The effective activation energy is plotted against the 
strain en — e. The solid line is a linear regression fit with 
slope 4.9 eV. 



IV. SIMULATION RESULTS 



The essential result of the MD simulations is the func- 
tional dependence of effective activation energy on lattice 
strain for diffusion on the GB plane. The function Q(e) 
is obtained from the Arrhenius plots of Diffusivity D|| 
vs. inverse temperature. This dependence is illustrated 
in Fig. 13 The solid line in Fig. Olis a linear regression fit 
to the data of the form Q(e) = 0.65eV - 4.9eVe. While 
experimental data with conditions similar to those sim- 
ulated here are not known to us, there are some general 
observations and comparisons that can be made with re- 
spect to the functional dependence of activation energy 
on strain. As expected, the zero-strain activation energy 
Q(0) — 0.65 eV is higher than that for most tilt bound- 
aries in copper, but significantly lower than the bulk ac- 
tivation energy of 1.3 eV for self-diffusion in copper. 

To place our data in the context of the analysis in Sec. 
Ill, let us compare the slope of the Q(e) data in Fig. |2| 
with the predicted slope EqV* (Eq. With En = 116 
GPa, we find V* = 0.60, where 17 = 1.2 x lO-^^m^ 
is the atomic volume for fcc-copper. While there is no 
experimental data with which to compare the value of 
V*, V*/il — 0.6 is consistent with vacancy-dominated 
diffusion and inconsistent with an interstitial-type mech- 
anism. This conclusion is based on the assumption that 
the migration volume contribution V*^ is negligible in 
comparison to the other volume terms in Eq. 1111 As- 
suming a vacancy mechanism in Eq. 1111 wc find that 
the relaxation volume satisfies V^/il — —1-2, while the 
assumption of a pure interstitial mechanism results in 

/ft — +4.8. The latter value is implausibly large, and 
we conclude that simple vacancies represent the domi- 
nant point defects in the GB. 
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V. DISCUSSION 

To estimate the technological significance of strain in 
interconnect device applications, consider that tensile 
stresses in copper interconnects are typically on the order 
of several hundred MPa, and the resulting strains depend 
on the crystallographic direction of the applied stress, but 
are on the order of 10""^. A reduction in effective acti- 
vation energy of « 5 eV per unit strain suggests that 
the strain effect is of order 10^^ eV, effectively doubling 
diffusivity at room temperature. While experimental dif- 
fusivity measurements can be uncertain to within orders 
of magnitude, static intrinsic strains typical of modern in- 
terconnect materials may significantly enhance diffusion 
near GBs in these materials. 

Any attempt to calculate activation energies from dif- 
fusion constants obtained with MD simulations is risky, 



in that the simulation time-scale may be shorter than 
the time taken for the defect concentration to reach 
equilibrium. However, the zero-strain activation energy 
Qo — 0.65 eV reported in this paper is within the range 
of values obtained by others for diffusion on a variety of 
GBs using molecular statics and assumptions about the 
dominant defect mechanism.t The evidence for a domi- 
nant vacancy-type mechanism is also consistent with sev- 
eral previous studies of defect formation on coherent twin 
GBs. 

Future computational work on grain boundary diffu- 
sion in copper will focus on measurements of V* for other 
technologically important GBs such as the S7 GB fam- 
ily. Like the S3 twin, E7 GBs are found in relatively high 
numbers in as-deposited copper thin films, but typically 
have significantly higher defect mobility than the S3 GB. 
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